Condensed Matter Physics, 2012, Vol. 15, No 2, 23001: 1-fTo] SMfflEHfflgHjD 
DOI: 10.5488/CMP.15.23001 JiJj^XWJBJi 
http://www.icmp.lviv.ua/journal 



Ion clusters and ion-water potentials 
in MD-simulationsI* 

Ph.A. BoppR K. Ibuklf 

1 Universitede Bordeaux, Department of Chemistry, 

351 Cours de la Liberation, Bat. A12, FR-33405 Talence cedex, France 

2 Doshisha University, Department of Molecular Chemistry and Biochemistry, Kyotanabe, Kyoto 610-0321, Japan 

Received December 31 , 201 1 , in final form March 1 3, 201 2 

A well known, if little documented, problem in many molecular simulations of aqueous ionic solutions at finite 
concentrations is that unrealistic cation-cation associations are frequently found. One might suspect a defect in 
the ion-ion interaction potentials, about which not much is known. However, we show that this phenomenon 
can also be traced to the fact that, in the pair-potential approximation, the cation-water potentials are too deep 
compared with the other ones and we investigate this phenomenon in some detail. We then attempt to draw 
some general conclusions. 
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1. Introduction 

It is a truism to stress the importance of aqueous salt solutions for many fields of science. Enormous 
efforts have consequently been made in the past decades to go beyond phenomenological considerations 
and to gain a fuller understanding of their structure and dynamics at the molecular level. A great many 
strategies in theoretical and computational physics and chemistry have been deployed for this purpose. 

It turns out that while the hydration of almost all simple positive ions has by now been investigated 
in Molecular Dynamics (MD) simulations at many levels of approximation (many pair potentials (see 
e.g. |1]), three body potentials |2], various QM-MM methods (see e.g. |3] and references therein), Car- 
Parinello MD (see e.g. |4]), Brownian dynamics and related approaches |5]), studies of solutions at finite 
concentrations are much scarcer. This is of course due to the much more difficult modeling task that one 
faces in this case. At least six interaction potentials (ten pair potentials in this approximation) must be 
determined in a consistent way. Their relative contributions to the total energy will vary widely when the 
salt concentration is varied. Any small imbalance between various terms can then lead to artefacts. We 
will come back to this point below. 

We have recently studied aqueous LiCl solutions in their entire concentration range at 300 K and 
normal densities ((J using a well established class of models that had been used for many salt solutions 
at low concentrations. Among the findings we observed at intermediate concentrations long-lived (with 
respect to the simulation times of a few 100 picoseconds) aggregations of Li + -ions. Figure[T]shows such a 
situation as it is found in a 10 m solution. Partly hydrated Li + -ions aggregate, leaving thus space for small, 
possibly interconnected, 'pools' of water with the weakly solvated Cl~-ions. The shape of these cationic 
aggregates is generally more or less linear and the ions are either in direct contact or barely solvent 
separated, as seen in the figure. 

The structural details of very concentrated ionic solutions at room temperature are little known. 
There is some evidence for the existence of, albeit possibly small (down to 5 H2O molecules), nano-pools 
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Figure 1. (Color online) Snapshot of a typical local arrangement in a 10 m LiCl solution at 300 K. Flexible 
water, a = 1, see text. The Li + -ions are dark green, the Cl _ -ions blue, water molecules closer than 2.5 A to 
a Li + -ion (hydration shell) are bright red, other water molecules pink. Hydrogen atoms are omitted. The 
Li + -clusters together with the hydration water and Cl~-ions close to Li + -ions (ion pairs) have been joined 
together by a green surface, the 'bulk' water by a brownish one. 

of water or worm-like structures in LiCl solutions at low temperatures |7]. Aggregations of some cations 
have also been found to influence the formation of certain complexes in aqueous solutions at room tem- 
perature d]. Even though they can thus not be entirely excluded, we do not know of any convincing 
experimental evidence for the existence of such aggregates, concerning either the ions or the solvent, in 
bulk solutions at higher temperatures. We therefore, decided to consider the Li + aggregates mentioned 
above as a possible artefact of our model and to study them in some detail in an attempt to trace their 
existence to some feature of our interaction model. We report here the results of this analysis and new 
simulations of the same solutions where this phenomenon has been eliminated by various means. 

2. Interaction models, general considerations 

The interaction models used in simulations of salt solutions are usually built up by combining a model 
which has been developed and tested for pure water (e.g. ST2, MCY, SPC, SPCE, TIP4P, BJH) with ion-water 
potentials developed in various ways, mostly by fitting under different constraints (e.g., the charge of 
the ion and geometry and partial charges of the water model) empirical functions to ab-initio energies. 
Simple charged Lennard-Jones (LJ) spheres are mostly, also in this work, assumed for the ion-ion interac- 
tions. 

Besides the problem discussed in the introductory section, several other ones have been identified. 
There are many critical discussions of the models for pure water [9-13], but we shall not enter this vast 
debate here. Pusztai et al. 11411 have in particular concluded from comparisons of scattering experiments, 
reverse Monte Carlo (rmc), and MD simulations that, at higher salt concentrations, the water-water inter- 
actions developed for pure water may no longer be adequate. 

There is indeed no guarantee whatsoever that interaction potentials of so different origins, often 
involving different approximations (effective potentials), can be combined in (almost) arbitrary propor- 
tions, depending on the salt concentration. Studying single ions, or solutions at low concentration, allevi- 
ates this problem. However, most 'real' solutions (biological fluids, sea water, in industry) are such that 
concentration effects cannot be neglected. 

In particular, an imbalance between the cation-water and water-water interaction potentials could 
lead to artefacts above a certain concentration, i.e. when the ion-water energies are no longer small 
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compared to e.g. the water-water ones. It can then happen that it becomes energetically more favorable 
to accept some ion-ion repulsion (especially between monovalent ions) and somehow allow these ions 
to associate so that they can share the remaining (scarce) water molecules to form a common 'hydration 
shell' for some sort of M" + aggregate. We explore this phenomenon (see e.g. in (§]) here by modifying the 
Li + -water interaction employed in |6] in two different ways: First by reducing the Coulomb-terms in the 
Li + -0 and Li + -H pair potential, and secondly by rigidifying in different ways the water molecule. 

3. Simulation details 

The details of the simulation are as in previous work @] except that either the Li + -water interac- 
tions are explicitly modified, or the geometry of the flexible water molecule is fixed, thus removing any 
mechanical polarization that may exist, see e.g. figure 12 in |6] for a distribution of the water molecular 
dipoles in LiCl solutions. Consequently, two options were pursued: a) to more or less arbitrarily reduce the 
Coulomb part of the Li + -oxygen and Li + -hydrogen pair potentials and b) to rigidify the water molecules 
in a reasonable geometry, see for comparison figure 7 in reference [15] and, as mentioned, figure 12 in ((J 
for the distributions of angles and the molecular dipole moments of flexible water in solutions. 
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Figure 2. (Color online) Li + -water energies as a function of the Li + -0 distance r, the parameter a (a- 
1.0,0.9,0.85,0.8 from bottom to top), and the geometry of the (flexible) water molecule (r H = 0.9572 A, 
^HOH = 104.5°, red curves; r H = 0.9721 A, Z H OH = 102.3°, green curves; / H = 0.9872 A, Z H OH = 101.7°, 
blue curves; rgn = 1-0115 A, Zhoh = 100.0°, pink curves, the curves have been shifted by 2 A in r- 
direction for better visibility.), for C2„ arrangements. The arrow shows that the potential at a = 0.9 for a 
deformed water molecule corresponds to the full potential for water in its approximate gas phase geom- 
etry. 

In all other respects, the simulations were identical to the previous ones: 556 BJH 1 16] water molecules 
and, depending on concentration, 10 to 150 ion pairs at experimental densities and 300 K. We have sim- 
ulated 1 m, 2 m, 4 m, 5 m, 8 m, 10 m, 12 m, and 15 m solutions. The simulations were run essentially 
under (TWFl-condifions, except that a Kast-type thermostat 1 17] was very loosely coupled (in order to 
minimize the perturbations of the dynamics, see below) to compensate an unavoidable numerical error 
during the s=10 6 integration steps. Ewald summation was used throughout. Data for the analysis were 
collected, after extensive equilibration runs, for 125 ps or 250 ps. The average temperature of the runs 
was (T) = (298+1.5) K. In one test case we increased the box size by a factor of 27 (15012 water molecules 
and 2700 ion pairs, 10 m solution). Starting with the end configuration of a smaller run we could extend 
this simulation only for about 6.75 ps. No changes in the structural data reported below could be detected 
in this run. 

The reduction in the cation-water interaction was simply achieved by multiplying the Coulomb terms 
of the Li + -0 and Li + -H pair potentials by factors a < 1, which is tantamount to multiplying for these 
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interactions the partial charges of Li + , and H by y/a. While this is inconsistent, e.g., from a dielectric 
standpoint (e.g.; which charges should be used to compute the dielectric constant e?), this is sufficient for 
our purposes here. We have applied factors of 0.8 a 1.0. Figure [2] shows the energies for a Li + -water 
supermolecule with C2„ geometry obtained for various values of a and for rigidified flexible BJH water 
molecules. We have also conducted simulations with such rigidified water (always with full charges, a = 
1) using two slightly different geometries: r 0H = 0.9572 A, Z H oh = 104.52° and Z H oh = 109.43°. 

4. Results and discussion 

Figure[3]shows for the flexible models some of the contributions to the total potential energies, divided 
by the salt concentration for better comparison, as a function of concentration, for the flexible water sim- 
ulations, with a as a parameter. The middle panel shows how lowering a, i.e. diminishing the Coulomb 
interaction between the Li + -ion and the water, changes this energy component in the simulations as a 
function of concentration. The top and bottom panels show how the average water-water and Cl~-water 
energies react to this modification of the Li + -water interaction. It is seen how lowering (increasing the 
magnitude) of the cation-water terms 'pushes up' the water-water energies, and that they can even be- 
come positive at higher concentrations for the larger as, as already noted previously (§)]. This means that, 
loosely speaking, the water molecules are, on average, unfavorably oriented with respect to each other 
(i.e. not adopting hydrogen-bond-type configurations), their orientations being presumably controlled by 
the cations. The CP-water energies (with unmodified pair potentials) increase in magnitude (i.e. are low- 
ered) when the Li + -water terms are diminished. Other energy terms are modified in a similar way, which 
will not be discussed in detail here. 

Figure H] shows, as an example, the influence of a on the water (O-O, O-H, and H-H radial pair dis- 
tribution functions (rdf)) for the 5 m solutions. The plot shows the rdf's and the 'running integration 
number' (number of neighbors) n{r) (both left ordinate). The modifications of these g-functions due to 
the modification of the Li + -water interaction are small, but clearly identifiable: The number of neighbors 
becomes slightly less when a goes from 1 to 0.8, although the effect seems not to be linear. 

To highlight the effect of the changes in the Li + -0 potential on all g-functions in the system, we have 
found it most instructive to plot the quantities 

n(r) 



Ara(r) = n rc (r)-rc(r) 

where rc and fc stand for 'reduced charges' (a=0.8, 0.85, 0.9) and 'full charges' (a = 1) in the Li + -water 
interactions, respectively. By analogy, rc will also be used for the rdf's obtained from rigid water simu- 
lations. Figure|4]also shows Arc(r) (right ordinate). The overall effect of a and concentration on An{r) is 
studied in a more systematic way in figures [6] 

Figure[5]is analogous to figure|4]for the Li + -0 and Li + -Li + functions. The effect of changing a is here 
much more pronounced and goes in the same direction as observed above for the water functions. It is 
seen in particular that the Li + -Li + aggregation is loosened: For a = 0.8, the number of next Li + -neighbors 
of a central Li + -ion is almost halved up to distances of about 8 A. 

Figure [6] shows An{r) as a function of r and the concentration for the three reduced charges. The 
panel (a) of this figure refers to the three g-functions of water, the panels (b) and (c) to the ion-water and 
the panel (d) to the ion-ion functions. Inspection of this figure and figures|4]and[5]reveals that modifying 
just the Li + -water interaction has a detectable effect on all radial pair distributions. This effect of a is 
generally most pronounced at intermediate concentrations (between about 6 and 12 m, say). As noted, 
the number of neighbors, i.e., the local particle density, is generally, at least up to intermediate distances, 
decreased when a goes from 1 to 0.8. The C1~-C1~ functions are an exception. We note here that the 
overall system densities were the same for all as ((A/V£)-simulations). 



= 47ipjg{r')r' 2 dr', 
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Figure 3. (Color online) Average water-water (top), 
Li + -water (middle), and Cl _ -water (bottom) ener- 
gies, divided by the molality of the LiCl-solution, as 
a function of the molality, for flexible water and 
values of a = 0.8, 0.85, 0.9, and 1.0. The arrows 
point in the direction of increasing a. 



Figure 4. (Color online) Water rdf's (go-o (top), 
gO-H (middle), and gn-H (bottom)) and their run- 
ning integration numbers n(r) (left abscissa, the 
labels give the values of g(r) and 0.5ra(r)) ; the 
variations of n{r) (An(r), right abscissa, see text) 
for the 5 m LiCI solutions. Flexible water, different 
a-values. 



Figure[7]shows that rigidifying the water molecules, i.e., reducing their average intermolecular ener- 
gies (due to the absence of the mechanical polarization of the flexible model) leads to similar, but smaller 
changes in the Li + -Li + rdf as explicitly reducing the Li + -water interaction. Similar observations are made 
for the other g-functions. 

Figure|8]shows the spectral densities of motion d{y) for the water molecule oxygens and the Li + -ions 
in the solutions of lowest and highest concentrations for the four values of a. The spectral density d{v) is 
the Fourier-cosine transform of the normalized velocity autocorrelation function c vv {t): 

I M N 

c vv {t) = c vv {t)/c vv {0) with C„„(£) = — — Y Y [v;(T/)-V,-(T/+ f)l, 

MN jri f-{ 
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f corr 

d(v)oc J c vv {t) cos (2nvt) dt, 
o 

where N is the number of equivalent particles (oxygen atoms, Li + -ions), M is the number of time origins 
(tj) of the correlation functions, v is the particle velocity and t con is the length of the correlation functions 
(not shown here), here generally 1 ps. Tests with other correlations lengths (i.e. different numbers M of 
correlation origins) yield the same results as the ones shown in the figure. 

Figure[8]shows that the effect of a on the translational dynamics is generally minor here; as expected, 
the spectral densities are (very slightly) shifted up for larger a values. The effect of the concentration, on 
the other hand, is considerable: The spectral densities of both species show well separated peaks at low 
concentration which merge more and more as the concentration increases. The v=0 term, which is pro- 
portional to the self-diffusion coefficient, also decreases with increasing concentration. This transition is 
more or less continuous when the concentration increases. For this reason, the densities for intermediate 
concentrations are not shown here. 



5. Further considerations 

The fact that the polarizabilities of the ion and of the water are not, or only in an average fashion and 
not explicitly, taken into account in the interaction models has often been criticized, especially if such a 
model is to be used in inhomogeneous environments, e.g. at surfaces |l£l|l9j]. However, another effect 
may also be important and of comparable, if not larger, magnitude: A minimal amount of charge transfer 
S\e\ (of the order of less than =;0.03 \e\ per water molecule in the first hydration shell of the central 
ion, say) will result in a discharged ion (in our example Li (1_6l5)+ , on the average) interacting with water 
molecules charged by the additional 5, which should be apportioned in some way to the partial charges of 
oxygen (-0.65971 e\ in the BJH model) and hydrogen (0.329851 e|). Of course, it is a priori not clear to which 
site of the model water the charge is best transferred. As figure [9] shows, such a transfer leads in a few 
plausible cases (e.g. S\e\ transferred entirely to the oxygen site of the water molecule, or 8\e\ distributed 
equally among the three water sites) to markedly shallower ion-water potential wells, comparable to 
what is seen in figure[2]for different values of a. Note that this plot does not take into account the fact 
that if several hydration water molecules are present, the ion will be even more discharged, as described 
above. 

In contrast to the present approach, the water-water interactions between first shell water molecules 
and the other water molecules would also be altered in case of a transfer of charge. The consequences 
of this will be explored in future work. In any case, these simple considerations show how effective 
a transfer of charge between an ion and its first shell of solvent molecules might be in modifying the 
interactions. 



23001-6 



Ion clusters and ion-water potentials 





r / A concentration (molality) r / A concentration (molality) 




Figure 6. (Color online) Variation of the running integration number as a function of r and salt concentra- 
tion, (An(r, concentration) , see text), for flexible water simulations and different values of a. The grey 
plane is the reference: a = 1, where A«(r) = 0. Panel (a): 0-0, 0-H, and H-H functions (water-water), 
panels (b), (c): Li + -0, Li + -H, Cl~-0, and Cl~-H functions (ion-water), panel (d): Li + -Li + , Li + -Cl~, and 
C1~-C1~ functions (ion-ion). Note the different scales for the ordinate in different rows. 



6. Summary 

We have studied aqueous LiCl solutions in the entire solubility range of this salt at room temperature 
and the densities corresponding to ambient pressure. Such systems are particularly interesting in terms 
of evaluating the underlying models because varying the concentration in a wide range means varying 
the contributions originating (in the pair potential approach) from various kinds of particle pairs (water- 
water, ion-water, ion-ion) in an equally wide range. It was the purpose of this work to show how small 
modifications of these interactions affect the results of the simulation. We chose first to vary the cation- 
water interaction since this is where we suspected in earlier work a bias leading to an artefact. Besides, 
we fixed the geometry of the flexible water model, thus eliminating polarization effects, which occur 
mainly in the ionic hydration shell, but also elsewhere. 

We have shown that decreasing the Li + -water infraction, in an admittedly arbitrary fashion, does 
indeed lead to qualitatively different results concerning the structure of the solutions, in particular at 
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Figure 7. (Color online) Selected panels like in figure [6] but for simulations with rigid water molecules: 
r 0H = 0.9572 A, Z H0H = 104.52° (left) and Z H0H = 109.43° (right). 



intermediate concentrations. Rigidifying the water, which leads to a reduction both in the ion-water and 
(less) the water-water energies, leads to similar, but smaller effects. 

This shows how sensitive the results of such simulations can be to details of the interaction model 
and to the balance between the components of such a 'global' model that may have been taken from very 
different sources and incorporating different approximations: effective potentials and ab-initio poten- 
tials, rigid vs. flexible models, polarization, charge transfer, etc.. The theoretical chemists and solution 
chemists have their work cut out. 




v/cm v/cm 

Figure 8. (Color online) Spectral densities of motion (Fourier-cosine transforms of the velocity autocorre- 
lation function) for the oxygen atoms of the water molecules (left) and the Li + -ions (right) in the 1 molal 
and 15 molal solutions for flexible water and the 4 different values of a. 
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Figure 9. (Color online) Similar to figure [2] Lithium-water potential wells for an assumed transfer of 
0.0125 |e| (green) and 0.0025|e| (blue) from a Li + -ion to a neighboring (hydration) water molecule. The 
fixed geometry of the (flexible) water molecule is tqh = 0.9572 A, Zhoh = 104.5°, as in the left curves in 
figure[2] and the red curve here is the one for this geometry and a = in figure[2] The lower (thin) green 
and blue curves here result form transferring all the charge from the ion to the oxygen atom of the water 
molecule, the upper curves (fat) when one third of the charge is transferred to each oxygen and hydrogen 
site of the water model. 
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IoHHi Ksiacrepi/i i noTemjia/ii/i ioH-BOfla b Mf\ ci/iMy/ummx 

OA Bond 31 , K. I6yi<P 

YHiBepcuTeT EopflO, Biflfli/i xiMii', FR-33405 Ta/iaHC, <J)paHL(in 

YHiBepcuTeT flOLuiiua, B\pj\\n MO/ieKy/wpHoT xi Mil' Ta 6ioxiMiT, KioTO 610-0321, flnoHifl 

flo6pe BiflOMOK), xona i Ma/io onwcaHOK), npo6/ieMOK) npi/i Bi/iKOHaHHi 4wce/ibHnx MO/ieKy/iapi-inx ci^My/ini4ii?i 
boahwx iOHHUx po34MHiB npi/i CKiHHeHHUx KOi-meHTpaLimx £ nonac™ Hepea^iCTUHHi acoi4iai4iT ™ny KaTiOH- 
KaTiOH. Mo>KHa npi/inyc™™, mo npunnHOK) i_(boro e flecfieKT iOH-iOHHi/ix noTei-mia/iiB B3aeMOfliT, npo flKuPi Bi- 
flOMO flOCTaTHbO Ma/io. Mi/i noKa3yeMO, luo u,e »BV\ui,e Moxe 6y™ cnpunnHeHe ™m, mo npn BMKopucTaHHi 
Ha6^wxeHHfl napHoT B3aeMOfliT noTei-mia/in KaTiOH-BOfla e 3aHaflTO mh6okwmh nopiBHHHO i3 iHLUHMW. Mn p,e- 
Ta/ibHO floaiiflxyeMO L4ePi ecfieicr y Haujifi po6oTi Ta HaMaraeMOCb ccfiopMyBaTH neBHi 3ara/ibHi bwchobkh. 

IOi K) h o b i cnoBa: CHMyjwLiii', MOJieKynapHa pyiHaMiKa, BO/jHO-iOHHi po3 L itiHii / ioHHa acoijiaLiin 
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